#include <cmath>
#include <iostream>
#include <string>

/**
 * @brief Sets material properties for elastic beam and bar elements (Type 401)
 * 
 * This function initializes material properties for elastic beam and bar elements.
 * It performs validation checks on material parameters and adjusts properties
 * for eccentricity if specified.
 * 

 */
void SETSG_201(std::vector<_Real_>& CM, _Real_ RO, int N, bool& ERROR) {
    // External dependencies
    extern std::vector<_Real_> PROP;  // Material properties from database
    extern std::string MSGBUF;        // Message buffer for error reporting

    // Local constants
    const std::string MSG = "*** ERROR *** PART ID";

    // Step 1: Set basic material properties
    CM[0] = RO;  // Density at index 0 (Fortran uses 1-based indexing)

    // Copy first 12 properties from PROP to CM (indices 1-12 in Fortran)
    for (int i = 0; i < 12; i++) {
        CM[i + 1] = PROP[i];
    }

    CM[15] = 0.0;  // Initialize index 16 (15 in C++ 0-based) to 0

    // Extract commonly used properties
    _Real_ YM = CM[1];   // Young's modulus (index 2 in Fortran)
    _Real_ ANU = CM[2];  // Poisson's ratio (index 3 in Fortran)
    _Real_ SEC = CM[5];  // Cross section area (index 6 in Fortran)
    _Real_ SSI = CM[7];  // Shear moment of inertia (index 8 in Fortran)
    _Real_ TTI = CM[8];  // Torsional moment of inertia (index 9 in Fortran)
    _Real_ RRI = CM[9];  // Rotational moment of inertia (index 10 in Fortran)

    // Step 2: Validate material properties
    ERROR = false;

    // Check Young's modulus
    if (YM <= 0.0) {
        MSGBUF = "<<==" + MSG + std::to_string(N) + 
                 " NONPOSITIVE YOUNG MODULUS =" + std::to_string(YM);
        std::cerr << MSGBUF << std::endl;
        ERROR = true;
    } else {
        // Calculate wave speed if Young's modulus is valid
        CM[19] = std::sqrt(YM / RO);  // Index 20 in Fortran
    }

    // Check Poisson's ratio
    if (ANU < 0.0 || ANU >= 0.5) {
        MSGBUF = "<<==" + MSG + std::to_string(N) + 
                 " POISSON COEFFICIENT NU =" + std::to_string(ANU);
        std::cerr << MSGBUF << std::endl;
        ERROR = true;
    }

    // Check cross section area
    if (SEC <= 0.0) {
        MSGBUF = "<<==" + MSG + std::to_string(N) + 
                 " NONPOSITIVE CROSS SECTION =" + std::to_string(SEC);
        std::cerr << MSGBUF << std::endl;
        ERROR = true;
    }

    // Check moment of inertia consistency
    int ICHK = 0;
    if (SSI != 0.0) ICHK += 1;
    if (TTI != 0.0) ICHK += 2;
    if (RRI != 0.0) ICHK += 4;

    // All should be zero (bar) or all non-zero (beam)
    if (ICHK != 0 && ICHK != 7) {
        MSGBUF = "<<==" + MSG + std::to_string(N) + 
                 " INCONSISTENT BEAM/BAR FORMULATION\n"
                 "<<== MOMENTS OF INERTIA SHOULD BE ALL NONZERO FOR A BEAM\n"
                 "<<== MATERIAL OR ALL ZERO FOR A BAR MATERIAL";
        std::cerr << MSGBUF << std::endl;
        ERROR = true;
    }

    // Step 3: Handle eccentricity if specified
    _Real_ STI = PROP[35];  // Index 36 in Fortran
    CM[125] = STI;          // Index 126 in Fortran
    CM[126] = 0.0;          // Index 127 in Fortran
    CM[127] = 0.0;          // Index 128 in Fortran

    _Real_ SBAR = PROP[36];  // Eccentricity in s-direction (index 37 in Fortran)
    _Real_ TBAR = PROP[37];  // Eccentricity in t-direction (index 38 in Fortran)

    if (SBAR != 0.0 || TBAR != 0.0) {
        // Adjust properties for eccentricity using parallel axis theorem
        CM[7] = SSI + SEC * TBAR * TBAR;  // Adjusted shear moment of inertia
        CM[8] = TTI + SEC * SBAR * SBAR;  // Adjusted torsional moment of inertia
        CM[9] = RRI + SEC * (SBAR * SBAR + TBAR * TBAR);  // Adjusted rotational moment
        CM[125] = STI + SEC * SBAR * TBAR;  // Adjusted product of inertia
        CM[126] = SEC * TBAR;  // First eccentricity term
        CM[127] = SEC * SBAR;  // Second eccentricity term
    }
}